Thermal conductivity of epitaxial graphene nanoribbons on SiC: effect of substrate 
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We study the effect of SiC substrate on thermal conductivity of epitaxial graphene nanoribbons 
(GNRs) using the nonequilibrium molecular dynamics method. We show that the substrate has 
strong interaction with single-layer GNRs during the thermal transport, which largely reduces the 
thermal conductivity. The thermal conductivity characteristics of suspended GNRs are well pre- 
served in the second GNR layers of bilayer GNR, which has a weak van der Waals interaction with 
the underlying structures. The out-of-plane phonon mode is found to play a critical role on the 
thermal conductivity variation of the second GNR layer induced by the underlying structures. 



PACS numbers: 

Graphene and Graphene nanoribbons (GNRs) are 
thought to be ideal materials for nanoelectronics due to 
their outstanding electronic and thermal properties]!- 
6j. In the production, the graphene nanomaterials 
can be either prepared by mechanical exfoliation from 
graphite B, Hi or by epitaxial growth on SiC substrate 0, 
flQI ] . However, the mechanical method is quite delicate 
and time consuming, which makes it unapplicable in 
the industry. The epitaxial growth method is nowa- 
days commonly accepted to represent a viable method 
of controllable growth for the fabrication of high quality 
graphene wafers The electronic properties of epi- 
taxial grap hene on SiC substrate had been extensively 
studied [lOl. Il2l - ll4j . It was found the electronic properties 
of graphene can be well preserved in both the single-layer 
(SL) graphene on SiC (0001) (C-terminated) [l(| , and the 
second layer of bilayer (BL) graphene on SiC(0001)(Si- 
terminated) [12], [l3| , having great application potential of 
epitaxial graphene in the nanoelectronics. Since the heat 
removal is a crucial issue in the nanoelectronic industry, 
the thermal conduction property of epitaxial graphene 
and GNRs becomes particularly important to its appli- 
cation in the nanoelectronics. 

During the last two years, the thermal conductivity of 
exfoliated graphene on different substrates (Si02, Cu) 
has been extensively studied, where only weekly cou- 
pled graphene-substrate interaction exists [l5l-[l9j. Dif- 
ferent from the exfoliated graphene case, the epitax- 
ial graphene-substrate interaction is much more compli- 
cated, and t he g eometry can even be distorted by the 
substrate [lU, l20j ]. Thus the thermal conductivity of epi- 
taxial graphene and GNRs is expected to be very different 
from the exfoliated ones. 

In this work, we use the nonequilibrium molecular dy- 
namics (NEMD)[rj|2ll-|23} method to study the thermal 
conductivity of epitaxial GNRs on 4H-SiC (0001) and 
(0001) surface. On the (0001) surface, both covalently 
bonded and weakly coupled GNR-substrate interaction 



conditions that were observed by the experiments are 
considered. On the (0001) surface, we consider both 
the SL and BL GNR-substrate interaction cases. In the 
BL GNR, we concentrate on the thermal conductivity of 
the second layer since the first layer is expected to have 
similar thermal conductivity as the SL GNR. Two typi- 
cal GNRs, i.e., armchair GNR (AGNR) and zigzag GNR 
(ZGNR), are considered, and we refer to AGNR/ZGNR 
with N dimer lines in width as N-AGNR/N-ZGNR for 
convenient representation [24j . 

The 4H-SiC substrate is modeled with four alternating 
Si and C atomic layers. One Si-C layer at the bottom of 
the sample is fixed. GNRs are placed on top of the SiC 
substrate, with infinite length alone the X direction and 
finite width along Y direction. The smallest cell of the 
GNR-SiC system is of 10.12 nm in length (X direction) 
and 2.15 nm (around 1 nm) for the SiC (GNR) in width 
(Y direction), containing 2128 SiC atoms and 320 (368) 
4-ZGNR (8-AGNR) atoms, respectively. 

In the geometry optimization, periodic boundary con- 
ditions are applied both along the X and Y directions. 
We use the Tersoff[25j potential to describe the C-C and 
C-Si bonded interactions, and the non-bonded van der 
Waals interaction is described by the Lennard- Jones (L J) 
potential[26j], which is only nonzero after the Tersoff co- 
valent potential goes to zero. The coupling between the 
long-range LJ potential and the short-range Tersoff po- 
tential is described by a cubic spline functkm|27]. 

In the NEMD simulation, we employ the velocity Ver- 
let method to integrate equations of motion with a fixed 
time step of 1 fs. Fixed boundary condition is applied 
along the X direction, where the outmost layers of each 
end of GNR (SiC) are fixed. Next to the boundaries, the 
adjacent 1 nm-long GNR (SiC) layers are coupled to the 
Nose-Hoover (28) thermostats with temperatures 310 and 
290 K, respectively. The thermal conductivity of GNRs 
k is then calculated from the Fourier law, 
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where J is the heat flux from the thermostats to the sys- 
tem, which is obtained from the Green-Kubo relation[29j, 
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FIG. 1: Optimized structures of 8- AGNR on SiC substrate 
for 4 cases: A, single- layer GNR on SiC (0001); B, bilayer 
GNR on SiC (0001); C, single-layer GNR on SiC (0001) 
with weakly coupled interaction; D, single-layer GNR on SiC 
(0001) with covalently bonded interaction. Carbon atoms are 
represented by gray balls; Si, by yellow balls. 



l30j . VT is the temperature gradient in the length direc- 
tion, which is defined as VT = (T L ~ T R )/L, where T L , 
Tr are the temperature of thermostats at the two ends, 
and L is GNR length. S is the cross-section area. In this 
work, we choose the interplanar spacing of graphite 3.35 
A as the GNRs' thickness. Moreover, all results given in 
this paper are obtained by averaging about 5 ns after a 
sufficient long transient time (5 ns) when a nonequilib- 
rium stationary state is set up. 

The following 4 kinds of epitaxial GNRs on SiC are 
considered: A, SL GNR on SiC (0001); B, BL GNR on 
SiC(0001); C, SL GNR on SiC (0001) with the_weakly 
coupled interaction; D, SL GNR on SiC (0001) with 
the covalently bonded interaction. Before geometry op- 
timization, the initial distance between GNR and SiC 
surface was set to be 2.3 A in cases A and B, 3.0 A m 
case C, and 2.0 A in case D. In Fig. 1, we show the 
optimized structures of 8- AGNR on SiC for the above 4 
cases. As one can see, the SL GNRs can either covalently 
bonded or weakly coupled to the substrate. In the cova- 
lently bonded case (case A, D), the formed C-Si (C-C) 
covalent bond between GNR and substrate has a length 
around 2.15 (1.60) A in the middle region, being consis- 
tent with the bond length of the graphene-SiC system in 
previous reports [HI, EH HH . While, in the edge region 
the C-Si (C-C) bond length is only about 2.10 (1.50) A, 
shorter than the bond length in the middle region, indi- 
cating that the edge atoms make the interaction between 
GNR and substrate stronger than that between graphene 
and substrate. In the weakly coupled case (case C), the 
mean distance between GNR and SiC surface is 3.30 A , 
similar with the graphene case [32| . 

In the BL GNR (case B), while, the mean distance be- 
tween the first and the second layer is 3.51 A, obviously 
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FIG. 2: The height profile of GNRs on SiC along the length 
direction for cases A (a), B (b), C (c), and D (d). 



larger than the interlayer distance of the suspended BL 
graphene. The larger interlayer distance comes from the 
larger shear modules of GNRs[33( and the large ripples 
formed on first GNR layer due to its strong interaction 
with SiC surface [3~ij]. Compared with the BL graphene, 
the interlayer force in BL GNR is smaller due to its finite 
width, while its shear modules is much larger. Thus, un- 
like the BL graphene, the weak interlayer van der Walls 
force can't make the second GNR layer follows the large 
ripples of the first GNR layer, and the second GNR layer 
still keeps a relatively planar structure. The large inter- 
layer distance appears at the trough of the ripple in first 
GNR layer. 

The corresponding height profile of GNRs on SiC along 
the length direction is shown in Fig. 2. From the figure, 
the ripples of GNRs that covalently bonded with SiC are 
obviously larger than that of GNR weakly coupled with 
SiC. The smallest and largest ripples appear in case C 
and D, respectively, both of which are on SiC (0001). 
Moreover, in the BL GNR the ripple of the second GNR 
layer is obviously smaller than that of first GNR layer 
(Fig. 2(b)), which corresponds to a larger interlayer dis- 
tance as discussed above. 

Fig. 3 shows the width dependence of thermal conduc- 
tivity of both AGNRs and ZGNRs on SiC. As shown in 
the figure, the thermal conductivity of GNRs covalently 
bonded with SiC is very different from that weakly cou- 
pled with SiC. In the covalently bonded case, thermal 
conductivity of GNRs is very low (below 30 W/mK), 
and has little variation when the width gets larger than 
2 nm. This shows that the strong covalent Si-C/C-C 
bonds formed between GNR and SiC surface have de- 
stroyed the intrinsic thermal conductivity of GNRs. In 
addition, the thermal conductivity of GNR on SiC (0001) 
(case D) is distinctly smaller than that on the (0001) sur- 
face (case A). This is attributed to the larger GNR ripples 
formed on the (0001) surface than that on (0001) surface 
(Fig. 2), which would induce stronger phonon scatter- 
ing. Different from the covalently bonded cases, thermal 
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FIG. 3: Width dependence of thermal conductivity of both 
AGNRs (a) and ZGNRs (b) on SiC in cases A (blue uptri- 
angle), B (cyan righttriangle ), C (red circle), and D (black 
square), with comparison of the suspended GNRs (magenta 
diamond) . 



conductivity of GNRs weakly coupled with SiC is much 
higher (case C). However, we find that some C-Si cova- 
lent bonds are newly formed between the edge of SL GNR 
and SiC(OOOl) surface during the NEMD simulation (in- 
set of Fig. 3), which is attributed to the large atomic 
amplitude of GNR in the direction perpendicular to the 
SiC surface when the system is thermostated to 300 K. 
The covalent bonds would induce additional edge local- 
ized phonon scattering on the GNR, and thus reduce the 
thermal conductivity. Since the number of newly formed 
covalent bonds is proportional to the length of GNR, 
compared with that of the suspended GNR, we expect 
the thermal conductivity will be largely reduced when 
the GNR length gets to /j,m-scale (experimental length). 

From Fig. 3, it is definite that, the thermal conduc- 
tivity of the second GNR layer in a BL GNR (case B) 
is most close to that of the suspended GNR. Also, the 
width-dependence of thermal conductivity is very simi- 
lar to each other: With the width increasing, the thermal 
conductivity of AGNR monotonously increases, while the 
thermal conductivity of ZGNR increases first and then 
decreases @. This implies the weak van der Waals inter- 
action from underlying structures (both the first GNR 
layer and SiC substrate) does not break the intrinsic ther- 
mal conductivity of the second GNR layer. 

In Fig. 4 we show the length dependence of thermal 
conductivity of the second GNR layers of BL GNR with 
comparison of the suspended GNRs. Similar with that 
of the suspended GNRs and carbon nanotubes frl 13H] . the 
thermal conductivity of the second 8-AGNR (4-ZGNR) 
layer monotonously increases with the length increasing 
and follows a power law of k ~ L' 9 , with ,3=0.38 (0.49). 
The value of (3 is very close to that of the suspended 
GNR (/3=0.32 (0.50)), showing that the ballistic and ID 
thermal transport characters [35t437j | are well preserved in 
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FIG. 4: Thermal conductivity k vs the length L in log-log 
scale for the second 8-AGNR (a), 4-ZGNR (b) layers in BL 
GNR, with comparison of the suspended GNRs. In both 
cases, k ~ L 13 , with j3 being close to that of the suspended 
GNRs. 



TABLE I: Thermal conductivity « of 8-ZGNRs with/without 
the out-of-plane vibrational constraint. The GNR length is 
kept at 10 nm. 



GNR type 


Constraint 


k (W/mK) 


Suspended SL GNR 


Free 


81.1 


Suspended SL GNR 


Constraint 


77.7 


2nd layer of BL GNR 


Free 


71.3 


2nd layer of BL GNR 


Constraint 


79.1 



the second GNR layer. 

To further explore the substrate effect on thermal con- 
ductivity of the second GNR layer of BL GNR on SiC, 
we freeze the out-of-plane atomic vibration of the second 
GNR layer and recalculate the thermal conductivity with 
comparison of the suspended SL GNR. For simplicity, 
here we only present the calculated thermal conductivity 
of ZGNRs, similar result is also obtained for AGNRs. As 
shown in Table 1, freezing the out-of-plane atomic vibra- 
tion decreases the thermal conductivity of the suspended 
SL GNR, while it considerably increases the thermal con- 
ductivity of the second GNR layer of BL GNR. 

This illustrates two facts. One is that the out-of- 
plane phonon mode in suspended SL GNR has a posi- 
tive contribution to the thermal transport, although it 
is not dominating. The other is that the thermal con- 
ductivity reduction of the second GNR layer of BL GNR 
mainly comes from the coupling between its out-of-plane 
phonon mode and the phonon modes of the underlying 
structures (including phonon modes of both first GNR 
layer 38| and SiC surface). It should be mentioned that, 
the thermal conductivity of the second GNR layer (79.1 
W/mK) in BL GNR is even higher than that of the sus- 
pended SL GNR (77.7 W/mK) when the out-of-plane 
atomic vibration is frozen. This interesting phenomenon 
indicates an increment of thermal conductivity can be re- 
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alized through the substrate coupling without the out-of- 
plane phonon. Recently, through a coupled atomic chain 
model, we have clarified that there exists a competitive 
mechanism on thermal conductivity in a coupling system: 
the phonon resonance effect that decreases thermal con- 
ductivity and phonon-band-up-shift effect that increases 
thermal conductivity [23j . In this case, the phonon res- 
onance effect mainly comes from the coupling between 
out-of-plane phonon mode of the second GNR layer and 
the phonon modes of the first GNR layer and SiC sub- 
strate. When the out-of-plane atomic vibration is frozen, 
the phonon resonance effect would be largely reduced, 
and the phonon-band-up-shift effect become dominated. 
Thus the thermal conductivity can be increased by the 
coupling. The results further confirm the existence of 
two competitive effects between the thermal conductive 
material and substrate. 

In summary, we have investigated the thermal conduc- 
tivity of epitaxial SL and BL GNRs on SiC substrate 
using the NEMD method. For the SL GNR, both co- 
valently bonded and weakly coupled GNR-substrate in- 
teraction conditions that were observed by experiments 
are considered. The thermal conductivity of SL GNRs 



in the covalently bonded condition is particularly low, 
while it is much higher in the weakly coupled condition. 
However, there appear some C-Si covalent bonds between 
the edge of SL GNR and SiC surface in the weakly cou- 
pled condition during the thermal transport, which in- 
duces additional edge localized phonon scattering and 
is expected to largely reduce the thermal conductivity 
when the GNR length gets to fim-scale. The second 
GNR layer of BL GNRs is found to have the highest 
thermal conductivity among all the epitaxial GNRs, and 
keeps much the same thermal conductivity characteristics 
as the suspended GNR. We find that the out-of-plane 
phonon mode of the second GNR layer plays a critical 
role on the thermal conductivity variation induced by 
the underlying structures. The existence of two com- 
petitive effects between thermal conductive material and 
substrate is further confirmed. We expect the present 
study can be helpful to the forthcoming applications of 
epitaxial graphene nanomaterials in the nanoelectronics. 
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